//Ghilardi and Zilberman (2024 JPE Macro)

// Endogenous Variables

var C p Div DivTau W Rk K N I q Y TauK TauD TauI theta A phi YY CC II KK qq pp DivDiv DivTau_annual IK IY IY_a AA theta_annual TauD_annual phi_annual N_annual;

// Exogenous Variables (Shocks)

varexo epsA epstheta epsTauD ;

// Parameters and Calibration

parameters betta delta alp N_ss theta_ss TauD_ss TauN_ss TauI_ss TauK_ss rhoA rhoTauD rhotheta A_ss gamma sd_theta sd_A 
phi_ss q_ss K_ss I_ss Y_ss C_ss h_ss W_ss Rk_ss Div_ss DivTau_ss p_ss I_y_ss K_y_ss I_k_ss C_y_ss Div_y_ss DivYtau;
 
betta=0.96;
delta=0.083;
alp=0.30;
N_ss=0.3;
theta_ss=0.0914;
TauD_ss=0.16;
TauN_ss=0; 
TauI_ss=0.0549;
TauK_ss=0.35;
rhoA=0.9; 
rhoTauD=0; 
rhotheta=0;
A_ss=1;
gamma=0.54;
sd_theta=0.01;
sd_A=0.0165;

//Steady State Model

//STEADY STATE LAGRANGE MULTIPLIER ON BORROWING CONSTRAINT

phi_ss=(delta/theta_ss)-(1-TauD_ss)*(1+TauI_ss);       //Use when considering a binding steady-state
//phi_ss=0;                                            //Use when the credit constraint is always slack. 

//REST OF STEADY STATE VALUES

q_ss=phi_ss+(1-TauD_ss)*(1+TauI_ss);
K_ss=(((alp*(1-TauK_ss))/((1+TauI_ss+((phi_ss)/(1-TauD_ss)))*((1/betta)-1)+delta*(1+TauI_ss)))^(1/(1-alp)))*N_ss;
I_ss=delta*K_ss;
Y_ss=A_ss*(K_ss)^(alp)*(N_ss)^(1-alp);
C_ss=Y_ss-I_ss;
h_ss=(1/C_ss)*(1-TauN_ss)*(1-alp)*(Y_ss/N_ss);
W_ss=(1-alp)*A_ss*((K_ss/N_ss)^alp);
Rk_ss=alp*A_ss*(N_ss/K_ss)^(1-alp);
Div_ss=(1-TauK_ss)*(Y_ss-W_ss*N_ss)-(1+TauI_ss)*I_ss;
DivTau_ss=(1-TauD_ss)*Div_ss;
p_ss=(betta*(1-TauD_ss)*Div_ss)/(1-betta);
I_y_ss=I_ss/Y_ss;
K_y_ss=K_ss/Y_ss;
I_k_ss=I_ss/K_ss;
C_y_ss=C_ss/Y_ss;
Div_y_ss=Div_ss/Y_ss;
DivYtau=(1-TauD_ss)*Div_y_ss;


model;

//NON-BINDING MODEL - Use these 2 equations when simulating the permanently slack model with phi_ss=0

//phi=0;

//q=(1-TauD)*(1+TauI+gamma*((I/K(-1))-delta));


//OBC MODEL - Use these 2 equations when allowing for OBC and phi_ss>0.  

phi=q-(1-TauD)*(1+TauI_ss+gamma*((I/K(-1))-delta));

0=min(phi,theta_ss*q*K(-1)-I);


// REST OF THE MODEL 

C=(1/h_ss)*(1-TauN_ss)*W;

p=betta*(C/C(+1))*((1-TauD(1))*Div(+1)+p(+1));

Div=(1-TauK)*(Y-W*N)-(1+TauI)*I-0.5*K(-1)*gamma*((I/K(-1))-delta)^2;

DivTau=(1-TauD)*Div; 

W=(1-alp)*A*((K(-1)/N)^alp);

Rk=alp*A*(Y/K(-1));

K=(1-delta)*K(-1)+I;

Y=A*((K(-1))^alp)*(N^(1-alp));

Y=C+I;

q=betta*(C/C(1))*((1-TauD(1))*((1-TauK(1))*Rk(1)+0.5*gamma*((I(1)/K)^2-delta^2))+q(1)*(1-delta+theta_ss*phi(1)));


//TAXES AND SHOCKS

TauK=TauK_ss;

TauI=TauI_ss;

TauD=(TauD_ss)-epsTauD;                                        //Div Tax shock process for perfect foresight solution

theta=((theta_ss)^(1-rhotheta))*((theta(-1))^rhotheta)*exp(sd_theta*epstheta);

A=((A_ss)^(1-rhoA))*((A(-1))^rhoA)*exp(sd_A*epsA);


//Annualized Equations

YY=100*(Y/Y_ss)-100;

CC=100*(C/C_ss)-100;

II=100*(I/I_ss)-100;

KK=100*(K/K_ss)-100;

qq=100*(q/q_ss)-100;

pp=100*(p/p_ss)-100;


DivDiv=100*(Div/Div_ss)-100;

DivTau_annual=100*(DivTau/DivTau_ss)-100;

IK=I/K(-1);       //deviations in levels

IY=I/Y;

IY_a=100*(IY/I_y_ss)-100;

phi_annual=phi;     //deviations in levels

N_annual=100*(N/N_ss)-100;

AA=100*(A/A_ss)-100;

theta_annual=100*(theta/theta_ss)-100;

TauD_annual=100*TauD;   //deviations in levels

end;


//// STEADY STATES DERIVED IN THE PARAMETER BLOCK

steady_state_model;

C=C_ss;
p=p_ss;
Div=Div_ss;
DivTau=DivTau_ss;
W=W_ss;
Rk=Rk_ss;
K=K_ss;
N=N_ss;
I=I_ss;
q=q_ss;
Y=Y_ss;
TauK=TauK_ss;
A=A_ss;
TauD=TauD_ss;
TauI=TauI_ss;
theta=theta_ss;
phi=phi_ss;
IK=I_ss/K_ss;
IY=I_y_ss;
YY=0;
CC=0;
II=0;
KK=0;
qq=0;
pp=0;
DivDiv=0;
N_annual=0;
phi_annual=phi_ss;
TauD_annual=100*TauD_ss;
theta_annual=0;
AA=0;

end;

steady;

shocks;

var epstheta; stderr 0; //Financial Shock NOT used in this code
var epsA; stderr 0;     //TFP Shock NOT used in this code

var epsTauD; periods 1:8; values 0.07;      //used for perfect foresight solution to determine the duration and size of the tax shock.   


end;

check;

steady;

simul(periods=50, stack_solve_algo = 0);

//USE FOR PERFECT FORESIGHT MODEL

%rplot YY;
%rplot CC;
%rplot II;
%rplot KK;
%rplot N_annual;
%rplot qq;
%rplot pp; 
%rplot DivDiv;
%rplot phi_annual;
%rplot TauD_annual;






